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Abstract 

Aims. We seek to find a shapelet-based scheme for deconvolving galaxy images from the PSF which leads to unbiased shear measurements. 
Methods. Based on the analytic formulation of convolution in shapelet space, we construct a procedure to recover the unconvolved shapelet 
coefficients under the assumption that the PSF is perfectly known. Using specific simulations, we test this approach and compare it to other 
published approaches. 

Results. We show that convolution in shapelet space leads to a shapelet model of order n h max = n g max + n f max with n f max and n g max being the 
maximum orders of the intrinsic galaxy and the PSF models, respectively. Deconvolution is hence a transformation which maps a certain 
number of convolved coefficients onto a generally smaller number of deconvolved coefficients. By inferring the latter number from data, we 
construct the maximum-likelihood solution for this transformation and obtain unbiased shear estimates with a remarkable amount of noise 
reduction compared to established approaches. This finding is particularly valid for complicated PSF models and low S /N images, which 
renders our approach suitable for typical weak-lensing conditions. 
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1. Introduction 

Shapelets have been proposed as an orthonormal set of two- 
dimensional functions to quantify shapes of galaxy images 
(Refregier 2003). They have several convenient mathematical 
properties suggesting their use in measurements of weak grav- 
itational lensing (Refregier & Bacon 2003 ). Their application 
to data with low signal-to-noise ratio, however, is hampered 
by a number of problems. First, their enormous flexibility al- 
lows shapelets to represent random noise patterns well, which 
generates their tendency to escape into high orders and to over- 
fit noisy images. Second, this same property not only affects 
the order, but also the scale of the best-fitting shapelet model. 
The scales of the shapelets decomposing noisy images is thus 
likely too high. These drawbacks are inherent properties of the 
shapelet method and have to be addressed for an effective us- 
age in astronomical image processing. At first sight, it seems 
natural to limit the order of the shapelets from above to avoid 
overfitting, and to limit also the scale from above to prevent 
shapelets from creeping into the noise. On the other hand, con- 
volution with the PSF alters the shape and spatial extent of 
imaged objects, and thus leaves an imprint on the maximum 
shapelet order and the scale size. 

Guided by these considerations, we address the follow- 
ing question: How should the maximum order and the spa- 
tial scale of a shapelet decomposition be determined in cases 
where PSF convolution significantly modifies the object's ap- 



pearance? This question aims at applications of shapelets in 
weak-lensing measurements, where the situation is particularly 
delicate because typically very small and noisy images are con- 
volved with structured, incompletely known PSF kernels with 
scales similar to those of the images. How can deconvolution 
schemes be constructed in this case in order to find significant 
and unbiased shear estimates? 

We show in Sect. [2] and in the Appendix that mathematical 
sum rules exist for the shapelet orders and the squared spa- 
tial scales of original image, kernel, and convolved image. We 
then proceed in Sect. [3] to devise an algorithm respecting these 
sum rules as well as possible, which leads to a deconvolution 
scheme based on the (possibly weighed) pseudo-inverse of a 
rectangular rather than the inverse of a quadratic convolution 
matrix. In Sect. |4| we demonstrate by means of simulations 
with different signal-to-noise levels and PSF kernels that our 
algorithm does indeed perform very well, and leads in most re- 
alistic cases to substantially improved results compared to pre- 
viously proposed methods. Our conclusions are summarized in 
Sect.EI 

2. Convolution in shapelet space 

A two-dimensional function /(x), e.g. a galaxy image, is de- 
composed into a set of shapelet modes by projection, 



/i 



■-r 



d 2 xf(x)B n (x;a), 
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where n = (ni,n 2 ) is a two-dimensional index and a is called 
scale size. The two-dimensional shapelet basis function 

B n (x;a) = a~ l cp ni (a~ l xi) cp n2 (a~ l x 2 ), (2) 

is related to the one-dimensional Gauss-Hermite polynomial 



cp n (x) = [2 n 7i2n\]-2 H n (x) e" 



(3) 



with H n (x) being the Hermite polynomial of order n. From the 
coefficients one can then reconstruct a shapelet model 



/(X): 



^ f n B n (x; a). 



(4) 



The number of shapelet modes - often expressed in terms of 
the maximum shapelet order n max = max(n\,n 2 ) - and the 
scale size have to be determined by an optimization algorithm 
( [Massey & Refregier|2005||Melchior et al.|2007| ), which min- 
imizes the modulus of the residuals f - f, or from empirical 
relations based on other measures of the object like FWHM or 
major and minor axes (Cha ng et al.|2004[ [Kuijken 2006 ). 
A convolution 



fc(x) = (/*#)(: 



%J-CK 



d 2 x'f(x')g{x-x'). 



(5) 



can be performed analytically in shapelet space. According to 
Eq. ([T]), the functions /, g and h are represented by sets of 
shapelet states f n , g n and h n with scale sizes a, ft and y. 

|Refregier| ([3003 ) and |Refregier & Bacon| ( [2003] ) showed 
that the coefficients of the convolved image h(x) are given by 



m,l 



(6) 



The 



where P nm = 2i C n ,m,igi is called convolution matrix. 
value of C n m i(a,fi, y) can be computed analytically. 

However, there is no clear statement on the scale size y and, 
in particular, on the maximum order n h max of the convolved ob- 
ject h. In appendix[A|we proof that the so-called natural choice 
( |Refregier|2003] ) 

y 2 = a 2 +p 2 (7) 

is indeed the correct choice for y and that the maximum order 
of the convolved object is given by 



n h = n f 

""max "max 



(8) 



While this result gives the highest possible mode of the con- 
volved object which could contain power, it does not tell us 
whether it does indeed have power, as this depends primarily 
on the ratio of scales a/p entering P n?m . This is demonstrated 
in Fig. [I] where we show the result of a convolution of a func- 
tion which is given by a pure B4 mode with a kernel represented 
by a pure B 2 mode. From this it becomes obvious that in a wide 
region around a/fi ^ 1 power is transfered to all even modes 
up to n = 6 (odd modes vanish because of parity, cf. Eq. ( |A.5| ) 
and the following discussion). If either a » ft or ft » a, the 
highest order of the larger object is also the highest effective 
order of the convolved object. Thus, we can generalize Eq. ([8]), 



''max 



(+1) a » p (kernel negligible) 

+ n 8 max a^p 

(+1) a <sc p (kernel dominant) 



(9) 




100.00 



Figure 1. One-dimensional convolution h = B^(x\ a)*B 2 (x; P). 
We plot the modulus of h n with n even (all odd modes have 
vanishing power), normalized by \h n \. 

where the option (+1) is taken if required by parity. 

For most cases, in particular for weak gravitational lensing, 
PSF and object scales are comparable, which means that we 
must not neglect the power transfer to higher modes. 

3. How to deconvolve 

While the correct values for y and n h max are clarified now, it is 
not obvious how these results need to be used in deconvolv- 
ing real data. We first comment on possible ways to undo the 
convolution and then discuss our finding in the light of mea- 
surement noise. 



3.1. Deconvolution strategies 

As Refregier & Bacon (2003) have already discussed, there are 
two ways to deconvolve from the PSF in shapelet space: 
- Inversion of the convolution matrix: According to Eq. ([6]), 
one can solve for the unconvolved coefficients, 



./m 



(10) 



- Fitting with the convolved basis system (Kuijken 1999 
Massey & Refregier 2005), which modifies Eq. (HI) such 



that it directly minimizes the residuals of h(x) w.r.t. its 
shapelet model 

n max 

m = YjUYj p «,mB m (x; a). (11) 

n m 

The second method is generally applicable but slow be- 
cause the convolution has to be applied at each iteration step 
of the decomposition process. On the other hand, the first ap- 
proach reduces deconvolution to a single step after the shapelet 
decomposition and is therefore computationally more efficient. 

According to Eq. ([8]), P is not quadratic and thus not in- 
vertible as suggested by Eq. fTO] ). In order to cope with this, 
we need to replace the inverse P~ l by the pseudo-inverse P^ = 
(P T P)~ l P T such that the equation now reads 



fm 



^ ^m,n n n- 



(12) 
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What seems as a drawback at first glance is effectively ben- 
eficial. Conceptually, this is now the least- squares solution of 
Eq. ([6J, recovering the most probable unconvolved coefficients 
from the set of noisy convolved coefficients. The underlying as- 
sumption of Gaussian noise in the coefficients h n holds for the 
most usual case of background-dominated images for which 
the pixel noise is Gaussian. 

Correlations among the coefficients h n , which arise from 
non-constant pixel weights or pixel correlations, can be ac- 
counted for by introducing the coefficient covariance matrix 
£ = {(h n - (h n ))(h n - (h n )) T ), which alters Eq. ^ to read 

(13) 



Maximizing the likelihood for recovering the correct uncon- 
volved coefficients leads to the weighted pseudo-inverse 



Pi = (p T zpy l p T z, 



(14) 



which replaces P^ in Eq. fT2] ). 

However, both approaches (direct inversion, Eq. fTO] ), or 
least-squares solution, Eq. fl2] )) would fail if P was rank- 
deficient. Refregier & Bacon] ([2003 ) argued that convolution 
with the PSF amounts to a projection of high-order modes onto 
low-order modes and therefore P can become singular. This 
is true only for very simple kernels (e.g. the Gaussian-shaped 
mode of order 0) with rather large scales. In fact, Eq. ^ tells us 
that convolution carries power from all available modes of / to 
modes up to order r^ max > r/ max , hence P is generally not rank- 
deficient. In practice we did not have problems in constructing 
P~ l or P^ when using realistic kernels. We therefore see no hin- 
drance in employing the matrix-inversion scheme and will use 
it in the course of this paper. 



3.2. Measurement process and noise 

Up to here, we have discussed (de-)convolution entirely in 
shapelet space, where this problem is now completely solved. 
For the following line of reasoning, we will further assume 
that the kernel is perfectly known and can be described by a 
shapelet model. 

Critical issues still arise at the transition from pixel to 
shapelet space: There are no intrinsic values of d max and a, 
and even if they existed they would not directly be accessi- 
ble to a measurement. While the first statement stems from the 
fact that we try to model a highly complicated galaxy or stellar 
shape with a potentially completely inappropriate function set, 
the second statement arises from pixelation and measurement 
noise occurring in the detector. 

However, the pixelated version of the shape can be de- 
scribed by a shapelet model, with an accuracy which depends 
on the noise level and the pixel size. Consider for example a 
galaxy whose light distribution strictly follows a Sersic pro- 
file. Modeling the cusp and the wide tails of this profile with 
the shapelet basis functions would require an infinite number 
of modes. But pixelation effectively removes the central singu- 
larity of the Sersic profile and turns the continuous light dis- 
tribution into a finite number of light measures, such that it is 
in principle describable by a finite number of shapelet coeffi- 
cients. Pixel noise additionally limits the spatial region within 




Figure 2. Sketch of the effect of a convolution on the power of 
shapelet coefficients. The detailed shape of the curves is nei- 
ther realistic nor important, but typical shapelet models show 
decreasing coefficient power with increasing order n. As con- 
volution does not change the overall power of an object but dis- 
tributes it over more coefficients, the average S /N of shapelet 
coefficients is lowered. The noise regime (represented by the 
gray area) is constant in the case of uncorrected noise. 



which the tails of Sersic profile remain noticeable and hence 
the number of required shapelet modes. 

Consequently, shapelet implementations usually determine 
n max by some significance measure of the model (x 2 in Massey 
& Refregier|[2Q05l |Melchior et"aL]|2007] ) or - similarly 



fix 

at a value which seems reasonable to capture the general 
features of the shape (e.g. Refre gier & Ba con 2003; Kuijken 
[20061 ). 

Fig. [2] schematically highlights an important issue of a 
significance-based ansatz: When the power in a shapelet co- 
efficient is lower than the power of the noise, it is considered 
insignificant, and the shapelet series is truncated at this mode 
(in Fig. [2] f n may be limited to n < 2 and h n to n < 3). Since 
convolution with a flux-normalized kernel does not change the 
overall flux or - as the shapelet decomposition is linear - the 
total coefficient power but generally increases the number of 
modes, the signal-to-noise ratio S /N of each individual coeffi- 
cient is lowered on average. Thus, after convolution more co- 
efficients will be considered insignificant and therefore disre- 
garded. 

This is equivalent to the action of a convolution in pixel 
space, where some objects' flux is distributed over a larger area. 
If the noise is independent of the convolution, demanding a 
certain S /N threshold results in a smaller number of significant 
pixels. 

The main point here is that we try to measure h n from data 
and from this /„ by employing Eq. fT2[ ). But if we truncate h n 



too early - at an order n h max <?c n J max + n„ 



the resulting un- 



convolved coefficients /„ are expected to be biased even if the 
convolution kernel is perfectly described. The reason for this 
is that by truncating we assume that any higher-order coeffi- 
cient is zero on average while in reality it is non-zero, it is just 
smaller than the noise limit. Every estimator formed from these 
coefficients is thus likely biased itself. 
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In turn, if we knew n^ max , we could go to the order de- 
manded by Eq. ([8]) and the deconvolution would map many 
noise-dominated high-order coefficients back onto lower-order 
coefficients. This way, we would not cut off coefficient power 
and our coefficient set would remain unbiased. Unfortunately, 
this approach comes at a price: Firstly, the resulting shapelet 
models are often massively overfitted, and secondly, obtaining 
unbiased /„ requires the knowledge of d max . The first problem 
can be addressed by averaging over sufficiently many galax- 
ies, while the second one can indeed be achieved by checking 
the S/N of the recovered f n after deconvolution. The average 
number of significant deconvolved coefficients gives an indi- 
cation of the typical complexity of the imaged objects as they 
would be seen in a measurement without convolution but using 
the particular detector characterized by its pixel size and noise 
level. 



3.3. Unbiased deconvolution method 

The previous consideration guides us to set up a deconvolu- 
tion procedure which yields unbiased deconvolved coefficients. 
Again, we assume perfect knowledge of the kernel g in shapelet 
space. 

- Given the noise level and the pixel size of the images, we 
initially guess n f max 

- We set the lower bounds rft max > n g max + h j nax and y > p. 

- We decompose each galaxy by minimizing the decomposi- 



tion^ 2 under these constraints. A value of n h n 



> n°„ 



hj nax is used only if x 2 > 1 otherwise. This yields h n and y. 

- By inverting Eq. ([7]), we obtain a. 

- Using the maximum orders and scale sizes for /, g and h 
in addition to g n , we can form the convolution matrix P 
according to Eq. ([6]). 

- By forming P^ and applying Eq. ([12]), we reconstruct f n . 

- By propagating the coefficient errors from the decomposi- 
tion through the same set of steps, we investigate the num- 
ber of significant coefficients and should find nf max if our 
initial guess was correct. 

Given the demanded accuracy, it might be necessary to ad- 
just the guess h^ max and reiterate the steps above. For this ap- 
proach, it is inevitable to split the data set in magnitude bins as 
the best value for n? max clearly depends on the intrinsic bright- 
ness. Further splitting (according to apparent size or brightness 
profile etc.) may be advantageous, too. 



4. Deconvolution microbenchmark 

There exists a growing number of shapelet-based decomposi- 
tion and deconvolution approaches published in the literature. 
In this section we will show that the method proposed here is 
indeed capable of inferring unbiased unconvolved coefficients. 
Moreover, employing the least- squares solution given by Eq. 
([12]) results in a considerable noise reduction, which is to be 
expected from this ansatz. 

At first, we want to emphasize that the simulations we 
use in this section are highly simplistic. Their only purpose is 
to investigate how well a certain decomposition/deconvolution 




Figure 3. Example of simulated galaxies used in our testbed: 

(a) intrinsic galaxy model with a = 2 and flux equaling unity, 

(b) after applying a shear y = (0.1, 0), (c) after convolving with 
PSFb from Fig. [4] with J3 - 2. The bottom panels show (c) after 
addition of Gaussian noise of zero mean and variance (t\\ (d) 
moderate noise, cr n = 10" 4 , (e) high noise, cr n = 10" 3 . (f) is the 
shapelet reconstruction of (e). Colors have logarithmic scaling. 



scheme can recover the unconvolved coefficients. By under- 
standing the performance of different approaches, we acquire 
the knowledge for treating more realistic cases. 

4.1. The testbed 

The construction of simulated galaxy images is visualized in 
Fig. [3] As intrinsic function we use a polar shapelet model 
with /o,o = fz,o = c, where c is chosen such that the model 
has unit flux, a is varied between 1.5 and 4. Given its ring- 
shaped appearance, this model is not overly realistic but also 
not too simple, and circularly symmetric. We apply a mild 
shear of y = (0.1,0), thus populate coefficients of order < 4, 
and convolve with five different realistic kernels g (cf. Fig. [4]) 
in shapelet space (employing Eqs. ([7]) & ([8]) with 1.5 < ft < 6). 
The pixelated version of the convolved object is then subject to 
N realizations of Gaussian noise with constant variance. 

Each of these simulated galaxy images is decomposed into 



shapelets again, yielding h n , using the code by Melchior et al. 
P007| ), where the optimization is constrained by fixing either 



n max or T' or both. h n is then deconvolved from the kernel g. 



As a diagnostic for the correctness of the deconvolved 
coefficients, we estimate the gravitational shear from the 



quadrupole moments G*7 of the light distribution (Bartelmann 
|& SchneiderpOOT] ), 



7Q 



1 



2 X 2 



IG11-G22 + 21G12 



Gn + G22 



(15) 



where Go is computed as a linear combination of all available 
deconvolved coefficients (Berge 2 005 [ 



Melch ior et al.|2007| ). 



We investigate five different approaches which differ in the 
choice of rft max , r/ max or the reconstruction of a. The different 
choices are summarized in Tab. Q] 



Full is the method we propose here (cf. Sect. 3.3 ); for the 



following tests, we set n\ 



f f 
= nl 



= 4. Signific is a variant 



P. Melchior et al.: Deconvolution with Shapelets 



5 




® e 



Figure 4. The kern els used in our bench mark: (a) model of 
PSF2 from STEP1 ( [Heymans et"aE||2Q06| ) with n g max = 4, (b) 
model of PSF3 from STEP1 with n 8 max = 4; (c) Airy disk model 
with n 8 max = 6; (d) model from a ray tracing simulation of a 
space-bourne telescope's PSF with n 8 max = 8 and n 8 max = 12 
(shown here). Colors have logarithmic scaling. 

Table 1. Overview of the parameter choices of the investigated 
methods. r^ max is the order of the decomposed object, a the esti- 
mate on the intrinsic scale and h^ max an estimate on the intrinsic 
order of /. 



Name 


"'max 


a 




n f 

n max 


Full 


— Wrnax Mmax 


Vr 2 - 




h f 

n max 


Signific 


> n 8 

— n max 


Vr 2 - 




n f 

n max 


Same 


n 8 

11 max 


Vr 2 - 


P 2 


n 8 

"■max 


ConstScale 


"■max 


r 




n 8 

"max 


Nmax2 


2 


Vr 2 - 


1 


2 



of Full, which bounds the decomposition order by the kernel 
order because coefficients beyond that are often insignificant, 
but makes use of our guess on rJ max . 

Same is similar to the one used by Kuijken (2006) with two 
differences: As discussed above, we employ the matrix inver- 
sion scheme (Eq. fLO} since P is square for this method) instead 
of fitting the convolved shapelet basis functions, and in our im- 
plementation x 2 is minimized w.r.t. a continuous parameter y, 
while |Kuijken| ( |2006| finds the best-fitting y = 2 n ^J3 with some 
integer n. Without knowing the increase of shapelet orders due 
to convolution given by Eq. ([8]), this represents the best-defined 
deconvolution approach. 

Refregier & Bacon (2003) stated that the approach 
ConstScale delivers the best results in their analysis. Nmax2, 
however, is an approach inspired by the naive assumption that 
such a decomposition scheme catches the essential shear infor- 
mation without being affected by overfitting. 



4.2. Performance with moderate noise 

The first set of simulations comprise galaxy models with peak 
S IN between 45 and 220 with a median of « 90 (an example 
is shown in Fig. [3ji); for each value of a and ft we created N = 
100 noise realizations. These high S /N values are more typical 
for galaxy morphology studies rather than for weak lensing, but 
we can see the effect of the convolution best. In this regime, 
problems with the deconvolution method become immediately 
apparent. 

Considering Fig. [5] we can ascertain that Full, Signific 
and Same perform quite well while ConstScale and Nmax2 
are clearly in trouble. This is not too surprising: By construc- 
tion, Nmax2 truncates the shapelet series at n l } nax = 2 and hence 



Full 
Signific 

Same 

ConstScale 
Nmax2 




Figure 5. Recovered shear y and intrinsic scale size a in a 
moderate-noise simulation (cr n = 10" 4 , a = 3, PSFb) as func- 
tions of the kernel scale p. In each panel the horizontal dashed 
line shows the true value of the quantity and the vertical dotted 
line shows the true value for a as reference. Errorbars (which 
are often too small to be visible) exhibit the standard deviation 
of the mean of N = 100 realizations. For visualization pur- 
poses, each method is slightly offset along J3. Simulations with 
different PSF models or a are qualitatively equivalent. 



misses all information contained in higher-order coefficients. 
One has to recall that the sheared model already has d max - 4, 
after convolution with PSFb (n 8 max = 4) it arrives at = 8. 
Nmax2 tries to undo the deconvolution with less information 
than contained in both sheared model and kernel individually. 
This is an enormously underconstrained attempt and leads to 
unpredictable behavior. ConstScale assumes that a can be ap- 
proximated by y and hence a is almost a increasing function 
of f3 (see bottom panel of Fig. [5|. According to Eq. ([7]), this 
ansatz is only applicable if ft is negligible. For very small kernel 
scales, we can indeed see a tendency to converge to the correct 
solution, but for all other situations, this choice is manifestly 
non-optimal. Because of the clearly problematic behavior of 
ConstScale and Nmax2, we exclude these two methods from 
the further investigation. 

This situation is very similar for other choices of a and 
other PSF models. To work out the general trends of the three 
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Figure 6. Recovered shear y and intrinsic scale size a (in units 
of the true scale size a) in the moderate-noise simulations in 
dependence of the PSF models from Fig. [4] (subscripts denote 
n 8 max ). Each data point represents the mean of the quantity for 
all available values of a and j3 (in total 60 independent com- 
binations), errorbars show the standard deviation of the mean. 
For visualization purposes, each method is slightly offset hori- 
zontally w.r.t. the others. 



remaining methods, we average over all scales a and ft and plot 
the results in dependence of the PSF model. 

The top and middle panels of Fig. [^confirm that all remain- 
ing methods yield essentially unbiased estimates of the shear, 
although we notice a mild tendency of Same and Signific to 
underestimate y i . This indicates that truncation of the decom- 
position order n^ ax = n g max might be insufficient for high S /N 
images. The fact that this underestimation is absent at higher 
kernel orders confirms this interpretation. 

Within the errors, the recovered scale size a is rather unbi- 
ased (see the bottom panel of [6]). For Same and Signific, we can 
see a clear shift of a for PSFc. The reason for this lies in the 
large spatially extent and wide wings of the Airy disk model in 
combination with a low r^ max . Since the entries of P depend in a 
non-linear way on a, this shift affects the recovery of the shear 
and leads to slightly poorer results. 




PSFa 4 PSFb 4 PSFc 6 PSFd 8 PSFdi 2 

Figure 7. Analogous to Fig. ^ but for the high-noise simula- 
tions. 



From this initial simulation with moderate noise we can 
conclude that one should respect Eq. ([7]) and must not truncate 
the shapelet series of h n severely. 

4.3. Performance with high noise 

We now consider a realistic weak-lensing situation by increas- 
ing the noise level by a factor of 10, hence 4.5 < S/N < 22 (cf. 
Fig. [3^). To balance the increased noise, we also increase the 
number of realizations to N = 1000. 

Considering Fig.|7j we can confirm that also for very noisy 
images the shear estimates from these three methods are not 
significantly biased. However, for Same we can see a remark- 
able drop of the mean of y\ and a drastic increase of the noise 
in y\ and 72 with the kernel order. Both findings are probably 
related to the usage of P~ l instead of P^ when performing the 
deconvolution. In contrast to the two methods we are propos- 
ing here, Same uses fJ max - n 8 max (cf. Table h|. For the typical 
weak-lensing scenario - characterized by rc^ ax < n 8 max , where 
all methods create a substantial amount of overfitting, cf. Fig. 
[3f -, this assumes to find a higher number of significant de- 
convolved coefficients then are actually available. These addi- 
tional, noise-dominated coefficients impact on Q t j and yg (cf. 
Eq. ([T5])), therefore these quantities become rather noisy them- 
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PSFc 6 PSFd 8 



Figure 8. Distance in shapelet space R s between the mean de- 
convolved and the true intrinsic coefficients in dependence of 
the PSF model for the three methods Full (top panel), Signific 
(middle panel) and Same (bottom panel). The mean is com- 
puted by averaging over all available values of a and j3 (in to- 
tal 60 independent combinations). Shown are the results for the 
moderate-noise simulations (solid line) and the high-noise sim- 
ulations (dotted line). 

selves. Given the fact that those high-order coefficients contain 
mostly arbitrary pixel noise which does not have a preferred di- 
rection, they also tend to dilute the available shear information 
from the lower-order coefficients, which explains the drop in 
y i . The estimate for y 2 is not affected as its true value was zero 
anyway. 

The superior behavior of Full and Signific in these low S /N 
simulations can also be seen more directly. As measure of the 
decomposition quality, we calculate the distance in shapelet 
space between the mean deconvolved coefficients f n and the 
true input coefficients /„, 



R 2 s = Yj«fo-fJ 2 ' 



(16) 



Fig. [8] confirms that as long as the kernel order is small, all 
three method perform quite similarly. But when the kernel or- 
der increases, Same tries to recover a quadratically increasing 
number of deconvolved coefficients whose individual signifi- 
cance is lowered at the same time. On the other hand, Full and 




10 12 
S/N 

Figure 9. Recovered shear y and intrinsic scale size a (in units 
of the true scale size a) in the high-noise simulations in depen- 
dence of the S IN of the convolved galaxy. The binning (dotted 
lines) is defined by the octiles of the S /N distribution, therefore 
all bins contain the mean values of approx. 7 combinations of 
a, J3 for each PSF model, in total ^35 independent settings. 
The data are plotted at the center of the bins and the methods 
are slightly offset horizontally for visualization purposes. Color 
code is as explained in Fig. [6] 



Signific make use of the redundancy of the overdetermined co- 
efficient set, which is created by applying a rectangular matrix 
P in Eq. ([6]). As a direct consequence of computing the least- 
squares solution via P\ the higher the number of convolved 
coefficients and the lower the number of significant intrinsic 
coefficients, the better these intrinsic coefficients can be recov- 
ered from noisy measurements. This explains the decrease of 
R s with the kernel order for these two methods. 

However, Full does not perform perfectly as well. The bot- 
tom panel of Fig. [7] reveals a bias on a, independent of the 
PSF model. The reason for this is again overfitting. As Full 
goes to higher orders than Signific and Same, it is even more 
affected by the pixel noise. As the decomposition determines 
y by minimizing^ 2 , y tends to become larger because this al- 
lows the model to fit a larger (increasingly noise-dominated) 
area, which reduces the overall residuals and thus^ 2 . Signific 
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and Same behave similarly when the kernel order - and hence 
the decomposition order - becomes larger. 

In order to prevent the shapelet models from creeping into 
the noisy areas around the object, it seems useful to constrain y 
not only from below but also from above. In addition to a guess 
on ntnax, we therefore impose a constraint p < y < ^fi 2 + a 2 max . 
Inferring both should be feasible when investigating observa- 
tional data. 

As our simulations comprise galaxy models of varying S /N 
- the models for both / and g have unit flux, so the surface 
brightness of the convolved object h depends on a and p -, it 
is illustrative to present the deconvolution results in S /N bins. 
Fig. [9] confirms that the two methods we propose here are very 
robust against image degradation. This is remarkable as many 
weak-lensing pipelines (and also Same in this paper) suffer 
from an underestimation of the shear, which becomes increas- 
ingly prominent with decreasing S /N (Mas sey et al.|2Q07| . Our 
statement from above, in which we related this drop to the high 
number of insignificant coefficients obtained from a deconvo- 
lution using Same, is further supported by this figure. It is ob- 
vious that - independent of the kernel model - a low S /N in 
pixel space results in a low S /N in shapelet space. By obtain- 
ing the least-squares solution for the /„, Full and Signific boost 
the significance of the recovered coefficients and thus perform 
better in the low S /N regime. The reason why y\ from Full is 
consistently but insignificantly lower than the estimates from 
Same is still somewhat unclear. A possible reason is the gener- 
ally higher number of shapelet coefficients h n for Full and thus 
a more noticeable noise contamination. 



5. Conclusions 

Based on an analytic consideration of shapelet convolution, we 
have studied algorithms for the PSF deconvolution of galaxy 
images in shapelet space. The starting point are sum rules for 
shapelet convolution, showing that the intrinsic shapelet orders 
of PSF and image add in the convolved image, and that the 
squares of their scales are also added. We suggest an algo- 
rithm respecting these sum rules as well as possible in pres- 
ence of noise, whose central step is the deconvolution of the 
convolved image with the pseudo-inverse of the convolution 
matrix. Applications to simulated images have shown that our 
algorithm performs very well and in many cases noticeably bet- 
ter than previously suggested methods. We identify three main 
reasons for the improved performance: 

- As the sum rule for the shapelet order shows, the convolu- 
tion transports power to higher shapelet modes. The mean 
signal-to-noise ratio of the convolved coefficients is thus re- 
duced. Our reduction of the order during the deconvolution 
increases the signal-to-noise without the need of calibra- 
tion. 

- The sum rules typically require a much higher shapelet or- 
der than thex 2 minimization, in particular if the PSF model 
is structured. Many of the high-order coefficients are thus 
highly insignificant. The significance of the coefficients is 
re-established by the reduction to the order which is 
chosen such that the shapelet expansion contains only sig- 



nificant coefficients. Fortunately, n? max depends mainly on 
the signal-to-noise ratio of the galaxy, but only weakly: 
many different galaxy shapes can be deconvolved with the 
same h^ max . Thus, binning the galaxies into broad signal-to- 
noise bins will suffice. 
- The STEP-2 project (Mass ey et al. 2007| ) has shown 
that shear measurements generally depend strongly on the 
galaxy brightness. Our algorithm seems to have the ad- 
vantage of lowering the influence of pixel noise in a well- 
defined manner. 

We shall proceed to study the performance of our algorithm in 
shear measurements under realistic conditions. 
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Appendix A: Convolution scale and order 

For simplicity we restrict ourselves to the one-dimensional 
case. From Eqs. (p} & ^ it is apparent that 



YjfmSi f dx'B m {x'\a)Bi{x-x'',p). (A.1) 

m,l J 



k(x) : 



We define 7 m> /(jc; a,p) as the integral in Eq. \A.\) and decom- 
pose it into shapelets with scale size y and maximum order N, 



I mJ (x; a, p) = ^ c n B n (x\ y). 



(A.2) 



Considering Eqs. ([2]), ([3j & ( |A.1| ), we recognize that N can- 
not be infinite but is determined by the highest modes of the 
expansions of / and g, which we will call M and L, respec- 
tively. Restricting to these modes and dropping all unnecessary 
constants, we can proceed, 



r r x a 1 

I dx\x') M exp - — (jc - x'f exp 
E(-D L+1 (t)^'J dx'(x') M+i ™V 



(x - x') 



A 2 



2p 2 

(x - x') 2 (x') 2 



ip 2 



2a 2 
(A3) 



where we expanded (x - x') L in the last step. By employing Eq. 
([7]) and substituting x = x f - ^ x, we can split the exponential, 



x exp 



i=0 ^ ' 



X 

~2f 



exp 



r 2 ~ 2 

X 



2a 2 fi 



(A.4) 
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/ 2 \M+i 

Again, we expand yx + ^ x) , which yields the desired ex- 



pression 

Im, 



M+L-j 



exp 



2y 2 



' 2a 2 £ 2 



(A.5) 



. Apart from the 



where we inserted Cj = J dx x j exp 



omitted constants, the second line of |A.5| is the definition of 
B M+L _j(x\y) (cf. Eqs. §2§ & ([3])) which shows that the natural 
choice is well motivated. Moreover, as j runs from to M + i, 
we see that the maximum order N is indeed M + L, as we have 
claimed in Eq. ([8]). Since Cj = if j is odd, the only states with 
non-vanishing power have the same parity as M + L. 

In summary, the natural choice is inherited from the 
Gaussian weighting function in Eq. ^ and the maximum order 
the result of a product of polynomials. 
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